# ==============================================================================
# file name: analysis-effects.R
# authors: Bernhard Clemm
# description: produces 
# - main results figure
# - Bayes' factors
# - exploratory analysis
# - alternative compliance measure
# relies on 02-recoding-traces.R and 01-recoding-survey.R
# ==============================================================================

# SETUP ========================================================================

rm(list = ls())
basedir <- paste0(dirname(dirname(
  rstudioapi::getSourceEditorContext()$path)), "/")
codedir <- paste0(basedir, "code/")
datadir <- paste0(basedir, "data/")
figdir <- paste0(basedir, "output/figures/")

library(tidyverse)
library(broom)
library(ivpack)
library(BayesFactor)

data_wide <- read.csv(paste0(datadir, "processed/data_wide.csv"))

dvs <- c("att_str", "att_imp", "ft", "stupid",
         "under", "malvol", "perpol", "supp_compr", "ptc",
         "know", "knowfac", "well_feel", "well_act")

# FUNCTIONS ====================================================================

model_to_table = function(model){
  x = summary(model)$coefficients %>% 
    data.frame() %>%
    mutate(variable = row.names(.)) %>%
    setNames(c("est", "se", "t", "p", "coeff")) %>% 
    select(., c("coeff", "est", "se", "p")) %>% 
    mutate(lb = est - (1.96*se)) %>%
    mutate(ub = est + (1.96*se))
  
  return(x)
}

# MAIN =========================================================================

## Run regressions ####

models = list()

for (var in dvs[1:13]){
  
  print(var)
  for (cntry in c("PL", "US")){
    
    print(cntry)
    
    dat <- data_wide %>% filter(country == cntry)
    
      if ((grepl("well", var) == T) | (var == "know" & cntry == "US")){
        dv <- paste0(var, "_post")
        iv <- ""
      } else {
        dv <- paste0(var, "_post")
        iv <- paste0(" + ", var, "_pre")
      }
      
      for (j in 1:7){
        
        eq = c(
          paste0(dv, " ~ condition_num", iv), 
          paste0(dv, " ~ compl_meanabs", iv, " | condition_num", iv), 
          paste0(dv, " ~ compl_self_post", iv, " | condition_num", iv), 
          paste0(dv, " ~ condition_num + news_before_mean + news_before_mean*condition_num", iv), # traces prior news expo
          paste0(dv, " ~ condition_num + news_self + news_self*condition_num", iv), # self-reported prior news expo
          paste0(dv, " ~ condition_num + news_like + news_like*condition_num", iv), # traces prior news ideo
          paste0(dv, " ~ condition_num + news_like_self + news_like_self*condition_num", iv) # self-reported prior news ideo
        )[j] 
        
        print(eq)
        
        modelname = c("ITT Effects", 
                      "CATE Effects (Behavioral)",
                      "CATE Effects (Self-Report)",
                      "Prior News Exposure (Behavioral)", 
                      "Prior News Exposure (Self-Report)",
                      "Prior Like-Minded News Exposure (Behavioral)", 
                      "Prior Like-Minded News Exposure (Self-Report)")[j]
        
        models[[length(models) + 1]] <- ivreg(as.formula(eq), data = dat) %>%
          model_to_table() %>%
          mutate(model = paste0("Model ", j)) %>% 
          mutate(modelname = modelname) %>% 
          mutate(country = cntry) %>%
          mutate(var = var)
        
      }}
}

models = do.call(rbind, models)

## Apply FDR adjustment ####

models$p_adjusted <- p.adjust(models$p, method = "fdr")
models <- models %>%
  mutate(sig = as.factor(ifelse(p_adjusted < 0.05, 1, 0)))

## Figure ####

fig_dat <- models %>%
  filter((
    modelname == "ITT Effects" & grepl("condition", .$coeff)) |
      (modelname == "CATE Effects (Behavioral)" & coeff == "compl_meanabs") |
      (modelname == "CATE Effects (Self-Report)" & coeff == "compl_self_post") |
      (modelname == "Prior News Exposure (Behavioral)" & coeff %in% c(
        "condition_num:news_before_mean", "condition_num:news_before_mean")) |
      (modelname == "Prior News Exposure (Self-Report)" & coeff %in% c(
        "condition_num:news_self","condition_num:news_self")) |
      (modelname == "Prior Like-Minded News Exposure (Behavioral)" & coeff %in% c(
        "condition_num:news_like", "condition_num:news_like")) |
      (modelname == "Prior Like-Minded News Exposure (Self-Report)" & coeff %in% c(
        "condition_num:news_like_self","condition_num:news_like_self"))) %>% 
  mutate(
    coeff = dplyr::recode(
      coeff, "condition_num" = "Treatment", 
             "compl_meanabs" = "Change in exposure (behavioral)",
             "compl_self_post" = "Change in exposure (self-reported)",
             "condition_num:news_self" = "Exposure × Treatment", 
             "condition_num:news_like_self" = "Exposure × Treatment", 
             "condition_num:news_before_mean" = "Exposure × Treatment", 
             "condition_num:news_like" = "Exposure × Treatment", 
             "condition_num:news_self" = "Exposure × Treatment",
             "condition_num:news_like_self" = "Exposure × Treatment", 
             "condition_num:news_before_mean" = "Exposure × Treatment",
             "condition_num:news_like" =   "Exposure × Treatment")) %>% 
  mutate(varno = case_when(
    # ITT effects
    modelname == "ITT Effects" ~ 13, 
    # CATE effects
    modelname == "CATE Effects (Behavioral)" ~ 11, 
    modelname == "CATE Effects (Self-Report)" ~ 9, 
    # News exposure - self-report
    modelname == "Prior News Exposure (Self-Report)" & 
      coeff == "Exposure × Treatment" ~ 7, 
    # News exposure - behavioral 
    modelname == "Prior News Exposure (Behavioral)" & 
      coeff == "Exposure × Treatment" ~ 5, 
    # Like-minded news exposure - self-report 
    modelname == "Prior Like-Minded News Exposure (Self-Report)" & 
      coeff == "Exposure × Treatment" ~ 3, 
    # Like-minded news exposure - behavioral 
    modelname == "Prior Like-Minded News Exposure (Behavioral)" & 
      coeff == "Exposure × Treatment" ~ 1)) %>%
  mutate(country = dplyr::recode(
    country, "PL" = "More News", "US" = "No News")) %>% 
  mutate(var = dplyr::recode(
    var, 'know' = '(a) Self-perceived Knowledge', 
        'knowfac' = '(b) Actual Knowledge', 
        'ptc' = '(c) Participation', 
        'supp_compr' = '(d) Support for Compromise',
        'att_str' = '(e) Attitude Polarization:\nAttitude Strength', 
        'att_imp' = '(f) Attitude Polarization:\nAttitude Importance',
        'ft' = '(g) Affective Polarization:\nFeeling Thermometer', 
        'under' = '(h) Affective Polarization:\nLack of Understanding', 
        'stupid' = '(i) Affective Polarization:\nStupid', 
        'malvol' = '(j) Attribution of Malevolence',
        'perpol' = '(k) Perceived Polarization',
        'well_feel' = '(l) Well-being: Mental', 
        'well_act' = '(m) Well-being: Physical'))

dist = 0.25
s = 2
s2 = 0.3

fig_dat <- fig_dat %>%
  mutate(varno_r = case_when(country == "More News" ~ varno + dist,
                             country == "No News" ~ varno - dist)) 

ggplot() + 
  geom_vline(xintercept = 0, size = 0.3, 
             linetype = "solid", colour = "grey60") + 
  geom_segment(data = fig_dat, 
               aes(y = varno_r, yend = varno_r, 
                   xend = ub, x = lb, colour = country,
                   alpha = sig), 
               size = s, lineend = "round") + 
  geom_point(data = fig_dat, 
             aes(y = varno_r, x = est, 
                 fill = country,
                 alpha = sig),
             color = "white", size = s2) + 
  scale_y_continuous(
    limits = c(1,14),
    breaks = c(1:14),
    labels = c(
     "Treatment × Exposure",
      expression(bold("Model 4b: Congenial News (Behavioral)")),
      "Treatment × Exposure",
      expression(bold("Model 4a: Congenial News (Self-Report)")),
      "Treatment × Exposure",
      expression(bold("Model 3b: News (Behavioral)")),
      "Treatment × Exposure",
      expression(bold("Model 3a: News (Self-Report)")),
      "Change in exposure",
      expression(bold("Model 2b: CATE (Self-Report)")),
      "Change in exposure",
      expression(bold("Model 2a: CATE (Behavioral)")),
      "Treatment",
      expression(bold("Model 1: ITT")))
    ) +
  scale_x_continuous(limits = c(-50,50)) +
  scale_fill_manual(values = c("grey60", "goldenrod1")) + 
  scale_colour_manual(values = c("grey60", "goldenrod1")) + 
  scale_shape_manual(values = c(21, 8, 4)) + 
  scale_alpha_discrete(name = "", range = c(0.45, 1), guide = "none") +
  theme_minimal() + 
  facet_wrap(var ~ ., ncol = 3, ) +
  theme(legend.position = c(.6, .08), 
        legend.title = element_blank(),
        legend.text = element_text(size = 15),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(1, "lines"), 
        axis.title = element_blank(), 
        axis.text.y = element_text(size = 13),
        strip.text = element_text(size = 13))

ggsave(paste0(figdir, "effects.png"), 
       device = "png", width = 15, height = 20)

# BAYES' FACTORS ===============================================================

models_bayes = data_frame()
i <- 1

for (var in dvs[1:13]){
  
  print(var)
  
  for (cntry in c("PL", "US")){
    
    dat <- data_wide %>% filter(country == cntry)
    
    print(cntry)
    
    if ((grepl("well", var) == T) | (var == "know" & cntry == "US")){
      dv <- paste0(var, "_post")
      iv <- ""
      plus <- ""
    } else {
      dv <- paste0(var, "_post")
      iv <- paste0(var, "_pre")
      plus <- " + "
    }
    
    for (j in 1:5){
      
      eq1 <- c(
        paste0(dv, " ~ condition_num", plus, iv), 
        paste0(dv, " ~ condition_num + news_before_mean + condition_num*news_before_mean", plus, iv), 
        paste0(dv, " ~ condition_num + news_self + condition_num*news_self", plus, iv),
        paste0(dv, " ~ condition_num + news_like + condition_num*news_like", plus, iv), 
        paste0(dv, " ~ condition_num + news_like_self + condition_num*news_like_self", plus, iv) 
      )[j] 
      
      eq1_vars <- unlist(strsplit(eq1, split = " ~ | \\+ |\\*")) %>% unique()
      
      dat_clean <- dat %>%
        select(all_of(eq1_vars)) %>% na.omit
      
      modelname = c("ITT Effects", 
                    "Prior News Exposure (Behavioral)", 
                    "Prior News Exposure (Self-Report)",
                    "Prior Like-Minded News Exposure (Behavioral)", 
                    "Prior Like-Minded News Exposure (Self-Report)")[j]
      
      # this is from https://richarddmorey.github.io/BayesFactor/#regression
      
      if ((grepl("well", var) == T & j == 1) | (var == "know" & cntry == "US" & j == 1)) {
        lmBF1 <- generalTestBF(as.formula(eq1), data = dat_clean)
        lmBF01 <- 1/lmBF1
        lmBF01_bf <- as.data.frame(lmBF01)$bf
      } else {
        lmBF1 <- generalTestBF(as.formula(eq1), data = dat_clean)
        lmBF01 <- (lmBF1[length(lmBF1)-1] / lmBF1[length(lmBF1)])
        lmBF01_bf <- as.data.frame(lmBF01)$bf
      }
      
      models_bayes[i, "dv"] <- dv
      models_bayes[i, "country"] <- cntry
      models_bayes[i, "model"] <- modelname
      models_bayes[i, "BF01"] <- lmBF01_bf
      
      i <- i + 1
    }
  }
}

models_bayes <- models_bayes %>% 
  mutate(dv = dplyr::recode(
    dv,
    'know_post' = 'Self-perceived Knowledge', 
   'knowfac_post' = 'Actual Knowledge', 
   'ptc_post' = 'Participation', 
   'att_str_post' = 'Attitude Strength', 
   'att_imp_post' = 'Attitude Importance',
   'ft_post' = 'AP: Feeling Thermometer', 
   'under_post' = 'AP: Lack of Understanding', 
   'stupid_post' = 'AP: Stupid', 
   'malvol_post' = 'Attribution of Malevolence',
   'supp_compr_post' = 'Support for Compromise',
   'perpol_post' = 'Perceived Polarization',
   'well_feel_post' = 'Well-being: Mental', 
   'well_act_post' = 'Well-being: Physical')) %>%
  mutate(country = dplyr::recode(country,
                          "PL" = "More News", 
                          "US" = "No News"))

## Table ####

kable(models_bayes, 
      caption = "Bayes' Factors comparing models to baselines", 
      format = "latex", booktabs = T, escape = F, linesep = "",
      row.names = F, longtable = T, 
      col.names = c("Outcome","Study", "Model", "BF01")) %>%
  kable_styling(font_size = 10) %>%
  collapse_rows(columns = 1, valign = "middle", latex_hline = "major") %>%
  collapse_rows(columns = 2, valign = "middle", latex_hline = "linespace") %>%
  save_kable(., file = paste0(tabdir, "models_bayes.tex"))

## Descriptives ####

sum(models_bayes$BF01 > 3) / length(models_bayes$BF01)
sum(models_bayes$BF01 > 10) / length(models_bayes$BF01)
sum(models_bayes$BF01[models_bayes$model == "ITT Effects"] > 10) / 
  length(models_bayes$BF01[models_bayes$model == "ITT Effects"])

# EXPLORATORY ==================================================================

models = list()

for (var in dvs[1:13]){
  
  print(var)
  for (cntry in c("PL", "US")){
    
    print(cntry)
    
    dat <- data_wide %>% filter(country == cntry)
    
    if ((grepl("well", var) == T) | (var == "know" & cntry == "US")){
      dv <- paste0(var, "_post")
      iv <- ""
    } else {
      dv <- paste0(var, "_post")
      iv <- paste0(" + ", var, "_pre")
    }
    
    for (j in 1:3){
      
      eq = c(
        paste0(dv, " ~ condition_num", iv), 
        paste0(dv, " ~ condition_num + leftright + leftright*condition_num", iv), 
        paste0(dv, " ~ condition_num + edu_high + edu_high*condition_num", iv) 
     )[j] 
      
      print(eq)
      
      modelname = c("ITT Effects", 
                    "Prior Ideology", 
                    "Prior Education")[j]
      
      models[[length(models) + 1]] <- ivreg(as.formula(eq), data = dat) %>%
        model_to_table() %>%
        mutate(model = paste0("Model ", j)) %>% 
        mutate(modelname = modelname) %>% 
        mutate(country = cntry) %>%
        mutate(var = var)
      
    }}
}

models = do.call(rbind, models)

## Apply FDR adjustment ####
models$p_adjusted <- p.adjust(models$p, method = "fdr")
models <- models %>%
  mutate(sig = as.factor(ifelse(p_adjusted < 0.05, 1, 0)))

## Figure ####

fig_dat <- models %>%
  filter((modelname == "ITT Effects" & grepl("condition", .$coeff)) |
           (modelname == "Prior Ideology" & coeff %in% c("condition_num:leftrightright")) |
           (modelname == "Prior Education" & coeff %in% c("condition_num:edu_highLow"))) %>% 
  mutate(
    coeff = dplyr::recode(coeff, 
                          "condition_num" = "Treatment", 
                          "condition_num:leftrightright" = "Ideology (right-wing) × Treatment", 
                          "condition_num:edu_highLow" = "Education (low) × Treatment")) %>% 
  
  mutate(varno = case_when(
    modelname == "ITT Effects" ~ 5, 
    modelname == "Prior Ideology" & 
      coeff == "Ideology (right-wing) × Treatment" ~ 3, 
    modelname == "Prior Education" & 
      coeff == "Education (low) × Treatment" ~ 1
  )) %>%
  mutate(country = dplyr::recode(country,
                                 "PL" = "More News", 
                                 "US" = "No News"
  )) %>% 
  mutate(var = dplyr::recode(var,
                             'know' = '(i) Political Engagement:\nSelf-perceived Knowledge', 
                             'knowfac' = '(ii) Political Engagement:\nActual Knowledge', 
                             'ptc' = '(iii) Political Engagement:\nParticipation', 
                             'att_str' = '(iv) Attitude Polarization:\nAttitude Strength', 
                             'att_imp' = '(v) Attitude Polarization:\nAttitude Importance',
                             'ft' = '(vi) Affective Polarization:\nFeeling Thermometer', 
                             'under' = '(vii) Affective Polarization:\nLack of Understanding', 
                             'stupid' = '(viii) Affective Polarization:\nStupid', 
                             'malvol' = '(ix) Negative System Perceptions:\nAttribution of Malevolence',
                             'supp_compr' = '(x) Negative System Perceptions:\nNo Support for Compromise',
                             'perpol' = '(xi) Negative System Perceptions:\nPerceived Polarization',
                             'well_feel' = '(xii) Well-being: Mental', 
                             'well_act' = '(xiii) Well-being: Physical'))

dist = 0.25
s = 2
s2 = 0.3

fig_dat <- fig_dat %>%
  mutate(varno_r = case_when(country == "More News" ~ varno + dist,
                             country == "No News" ~ varno - dist)) 

ggplot() + 
  geom_vline(xintercept = 0, linetype = "dashed", colour = "grey60") + 
  geom_segment(data = fig_dat, 
               aes(y = varno_r, yend = varno_r, 
                   xend = ub, x = lb, colour = country,
                   alpha = sig), 
               size = s, lineend = "round") + 
  geom_point(data = fig_dat, 
             aes(y = varno_r, x = est, 
                 fill = country,
                 alpha = sig),
             color = "white", size = s2) + 
  scale_y_continuous(
    limits = c(0,6),
    breaks = c(1:6),
    labels = c(
      "Treatment × Education (low)",
      expression(bold("Model 5b: Heterogeneity by education")),
      "Treatment × Ideology (right-wing)",
      expression(bold("Model 5a: Heterogeneity by ideology")),
      "Treatment",
      expression(bold("Model 1: ITT")))
  ) +
  scale_x_continuous(limits = c(-50,50)) +
  scale_fill_manual(values = c("grey60", "goldenrod1")) + 
  scale_colour_manual(values = c("grey60", "goldenrod1")) + 
  scale_shape_manual(values = c(21, 8, 4)) + 
  scale_alpha_discrete(name = "", range = c(0.37, 1), guide = "none") +
  theme_minimal() + 
  facet_wrap(var ~ ., ncol = 3) +
  theme(legend.position = c(.6, .08), 
        legend.title = element_blank(),
        legend.text = element_text(size = 15),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(1, "lines"), 
        axis.title = element_blank(), 
        axis.text.y = element_text(size = 13),
        strip.text = element_text(size = 13))

ggsave(paste0(figdir, "effects_explore.png"), 
       device = "png", width = 15, height = 10)

# ALTERNATIVE COMPLIANCE =======================================================

## Run regressions ####

models = list()

for (var in dvs[1:13]){
  
  print(var)
  for (cntry in c("PL", "US")){
    
    print(cntry)
    
    dat <- data_wide %>% filter(country == cntry)
    
    if ((grepl("well", var) == T) | (var == "know" & cntry == "US")){
      dv <- paste0(var, "_post")
      iv <- ""
    } else {
      dv <- paste0(var, "_post")
      iv <- paste0(" + ", var, "_pre")
    }
    
    for (j in 1:3){
      
      eq = c(
        paste0(dv, " ~ condition_num", iv), 
        paste0(dv, " ~ pol_news_compl_meanabs", iv, " | condition_num", iv), 
        paste0(dv, " ~ pol_compl_meanabs", iv, " | condition_num", iv)
        )[j] 
      
      print(eq)
      
      modelname = c("ITT Effects", 
                    "CATE Effects (Behavioral: Political news)",
                    "CATE Effects (Behavioral: Political content)")[j]
      
      models[[length(models) + 1]] <- ivreg(as.formula(eq), data = dat) %>%
        model_to_table() %>%
        mutate(model = paste0("Model ", j)) %>% 
        mutate(modelname = modelname) %>% 
        mutate(country = cntry) %>%
        mutate(var = var)
      
    }}
}

models = do.call(rbind, models)

models <- models %>%
  mutate(sig = as.factor(ifelse(p < 0.05, 1, 0)))

## Figure ####

fig_dat <- models %>%
  filter((modelname == "ITT Effects" & grepl("condition", .$coeff)) |
           (modelname == "CATE Effects (Behavioral: Political news)" & coeff == "pol_news_compl_meanabs") |
           (modelname == "CATE Effects (Behavioral: Political content)" & coeff == "pol_compl_meanabs")) %>% 
  mutate(
    coeff = dplyr::recode(
      coeff, 
      "condition_num" = "Treatment", 
      "pol_news_compl_meanabs" = "Change in exposure (behavioral)",
      "pol_compl_meanabs" = "Change in exposure (self-reported)")) %>% 
  
  mutate(varno = case_when(
    
    # ITT effects
    modelname == "ITT Effects" ~ 5, 
    
    # CATE effects
    modelname == "CATE Effects (Behavioral: Political news)" ~ 3, 
    modelname == "CATE Effects (Behavioral: Political content)" ~ 1)) %>%
  mutate(country = dplyr::recode(country,
                                 "PL" = "More News", 
                                 "US" = "No News"
  )) %>% 
  mutate(var = dplyr::recode(var,
                             'know' = '(a) Self-perceived Knowledge', 
                             'knowfac' = '(b) Actual Knowledge', 
                             'ptc' = '(c) Participation', 
                             'supp_compr' = '(d) Support for Compromise',
                             'att_str' = '(e) Attitude Polarization:\nAttitude Strength', 
                             'att_imp' = '(f) Attitude Polarization:\nAttitude Importance',
                             'ft' = '(g) Affective Polarization:\nFeeling Thermometer', 
                             'under' = '(h) Affective Polarization:\nLack of Understanding', 
                             'stupid' = '(i) Affective Polarization:\nStupid', 
                             'malvol' = '(j) Attribution of Malevolence',
                             'perpol' = '(k) Perceived Polarization',
                             'well_feel' = '(l) Well-being: Mental', 
                             'well_act' = '(m) Well-being: Physical'))

dist = 0.25
s = 2
s2 = 0.3

fig_dat <- fig_dat %>%
  mutate(varno_r = case_when(country == "More News" ~ varno + dist,
                             country == "No News" ~ varno - dist)) 

ggplot() + 
  geom_vline(xintercept = 0, size = 0.3, 
             linetype = "solid", colour = "grey60") + 
  geom_segment(data = fig_dat, 
               aes(y = varno_r, yend = varno_r, 
                   xend = ub, x = lb, colour = country,
                   alpha = sig), 
               size = s, lineend = "round") + 
  geom_point(data = fig_dat, 
             aes(y = varno_r, x = est, 
                 fill = country,
                 alpha = sig),
             color = "white", size = s2) + 
  scale_y_continuous(
    limits = c(1,6),
    breaks = c(1:6),
    labels = c(
      "Change in exposure",
      expression(bold("Model 2a: CATE (Behavioral: Political news)")),
      "Change in exposure",
      expression(bold("Model 2a: CATE (Behavioral: Political content)")),
      "Treatment",
      expression(bold("Model 1: ITT")))
  ) +
  scale_x_continuous(limits = c(-50,50)) +
  scale_fill_manual(values = c("grey60", "goldenrod1")) + 
  scale_colour_manual(values = c("grey60", "goldenrod1")) + 
  scale_shape_manual(values = c(21, 8, 4)) + 
  scale_alpha_discrete(name = "", range = c(0.37, 1), guide = "none") +
  theme_minimal() + 
  facet_wrap(var ~ ., ncol = 3, ) +
  theme(legend.position = c(.6, .08), 
        legend.title = element_blank(),
        legend.text = element_text(size = 15),
        panel.grid.minor = element_blank(),
        panel.spacing = unit(1, "lines"), 
        axis.title = element_blank(), 
        axis.text.y = element_text(size = 13),
        strip.text = element_text(size = 13))

ggsave(paste0(figdir, "effects_compliance_alt.png"), 
       device = "png", width = 15, height = 10)



